# Replication code for "Supply-Side Reforms to Oil and Gas Production on Federal
# Lands: Modeling the Implications for CO2 Emissions, Federal Revenues, and Leakage"
# Author: Brian Prest

#****# Preamble ------------------------------------------------------------------------------------
"%&%" <- function(x,y) paste0(x,y)

library(data.table)
library(ggplot2)
library(lubridate)
library(tidyverse)
library(sp)
library(dplyr)
library(sf)
library(rnaturalearth)
library(car)
library(sandwich)
library(AER)

project.dir = '/path/to/replication code/'
output.dir = project.dir %&% '/Output/'
script.dir = project.dir%&%'/scripts/'
data.dir =   project.dir %&% '/Data/'
di.data.folder = data.dir %&% '/Drillinginfo/'
working.data.dir = data.dir %&%'/Working_Data/'

dir.create(output.dir, recursive = TRUE)
dir.create(working.data.dir, recursive = TRUE)

options(digits=7)

# RFF color palette
palette(c(rgb(4, 39, 60, maxColorValue=255), # RFF black
          rgb(116,100,94, maxColorValue=255), # RFF brown
          rgb(80,177,97, maxColorValue=255), # RFF green
          rgb(136,196,244, maxColorValue=255), # RFF blue
          rgb(255, 102, 99, maxColorValue=255), # RFF coral
          rgb(118, 94, 165, maxColorValue=255), # RFF purple
          rgb(244, 162, 95, maxColorValue=255), # RFF orange
          rgb(235, 211, 103, maxColorValue=255), # RFF yellow
          rgb(224, 228, 231, maxColorValue=255))) # RFF light gray

rep.row<-function(x,n) matrix(rep(x,each=n),nrow=n)
rep.col<-function(x,n) matrix(rep(x,each=n), ncol=n, byrow=TRUE)

# Read Data, Compute Dates and Well Types 
source(script.dir%&%'read_data_and_identify_wells.R')

# Make plots of historical production
source(script.dir%&%'plot_historical_production.R')

# Make map of well locations
source(script.dir%&%'map_wells.R')

# Make panel of spudcounts
source(script.dir%&%'create_spudcounts_panel.R')

# Estimate drilling elasticities
source(script.dir%&%'estimate_pooled_drilling_elasticities.R')

# Estimate spud-to-prod distributions
source(script.dir%&%'spud_to_prod_distributions.R')

# Calculate production profiles
source(script.dir%&%'prep_profiles.R')

# Set up price paths for WTI and HH
source(script.dir%&%'prep_prices.R')

# Prep parameters for simulation (shock date, base spuds, spud-to-prod times, production profiles)
source(script.dir%&%'prep_sim_parameters.R')

# Baseline production
source(script.dir%&%'project_well_level_production.R')

# Define simulation function
source(script.dir%&%'define_sim_functions.R')

# Helper functions. Plots simulation results. Not necessary.
source(script.dir%&%'sim_plotting_functions.R')

# Functions to extract results of simulations
source(script.dir%&%'simulation_summary_helpers.R')

### Run simulations ----
# Clear extraneous data, including full dataset (dt) and spudcount panels
dropvars = grep('spudcounts', ls(), v=T)
dropvars = dropvars[dropvars!='spudcounts.class']
rm(list=c('dt','ocean','over.out','locs','grp.idx', dropvars)); gc()

# Save working data for later use, without having to rerun all code above.
# Use this copy when running the alternative sensitivities below, such as changing
# the ROW supply elasticities.
save.image(working.data.dir%&%'/working_fed_lands_data_'%&%today()%&%'.RData')
# load(working.data.dir%&%'/working_fed_lands_data_2021-11-01.RData')

# These should all be the same for simulation to work properly
dimnames(base.spuds) # baseline spuds (monthly)
dimnames(lm.regs) # elasticity estimates
dimnames(stp.den)[-1] # spud-to-prod distributions
dimnames(roy.base)

### Set sensitivity parameters. Base case is:
#   dem.elasts.to.use = 'Base' # baseline demand elasticities
#   prices.to.use = 'futures' # base price forecast on futures strip
#   global.gas.market = TRUE # model gas market as global
#   low.supp.elast = FALSE # don't use low ROW supply elasticity sensitivity
#   us.supp.elast = FALSE # don't use high (US-based) ROW supply elasticity sensitivity

## Demand elasticities
dem.elasts.to.use = 'Base'
# dem.elasts.to.use = 'High'

## Price baseline
prices.to.use = 'futures'
# prices.to.use = 'IEA'

# Global gas market?
global.gas.market = TRUE
# global.gas.market = FALSE

## Low ROW supply elasticity sensitivity?
low.supp.elast = FALSE
# low.supp.elast = TRUE
if (low.supp.elast) {
  weo.projections[, row_crude_elasticity := 0.5*row_crude_elasticity]
  weo.projections[, row_gas_elasticity := 0.5*row_gas_elasticity]
}
## High supply elasticity sensitivity? Must run this after the baseline run has generated the 'US_avg_elasticities_futures_Base_ROW.csv' file
us.supp.elast = FALSE
# us.supp.elast = TRUE
if (us.supp.elast) {
  us.elasts = fread(output.dir%&%'/US_avg_elasticities_futures_Base_ROW-elasts-IEA.csv')
  weo.projections[, row_crude_elasticity:=NULL]
  weo.projections[, row_gas_elasticity:=NULL]
  setnames(us.elasts, c('date','row_crude_elasticity','row_gas_elasticity'))
  us.elasts[, date:=as.Date(date)]
  weo.projections = merge(weo.projections, us.elasts, by='date')
}

# Run simulations
for (global.gas.market in c(TRUE, FALSE)) { 
  for (prices.to.use in c('futures','IEA')) {
    for (dem.elasts.to.use in c('Base','High')) {
      if (us.supp.elast & !global.gas.market) next # don't need to run both sensitivities at the same time
      if ((us.supp.elast | !global.gas.market) & prices.to.use!='futures') next # only run the sensitivities with baseline prices
      if (dem.elasts.to.use=="High" & global.gas.market==FALSE) next 
      
      message(paste('Iteration start time:', Sys.time()))
      
      source(script.dir%&%'run_sims.R')
      
      if (prices.to.use=='futures' & dem.elasts.to.use=='Base') file.nm.append = '_'%&%today()%&%'.RData'
      if (prices.to.use=='IEA' & dem.elasts.to.use=='Base') file.nm.append = '_high_prices_'%&%today()%&%'.RData'
      if (prices.to.use=='futures' & dem.elasts.to.use=='High') file.nm.append = '_high_elasts_'%&%today()%&%'.RData'
      if (prices.to.use=='IEA' & dem.elasts.to.use=='High') file.nm.append = '_high_elasts_and_prices_'%&%today()%&%'.RData'
      if (us.supp.elast) file.nm.append = '_HighUS_ROW_elasts_'%&%today()%&%'.RData'
      if (!global.gas.market) file.nm.append = '_Domestic_gas_market_'%&%today()%&%'.RData'
      save.image(working.data.dir%&%'/working_fed_lands_data_modelrun'%&%file.nm.append)
      message(paste(prices.to.use%&%' price', dem.elasts.to.use%&%' demand elasts', 'ROW High Elast Sens? '%&% us.supp.elast, 'Global gas market? '%&%global.gas.market, 'done '%&% Sys.time(), sep=', '))
      
    }
  }
}
working.data.loc = working.data.dir%&%'/working_fed_lands_data_modelrun_2021-11-01.RData'
working.data.loc.high.e = working.data.dir%&%'/working_fed_lands_data_modelrun_high_elasts_2021-11-01.RData'

# Can load saved copies of previous run models under alternative assumptions, as follows,
# and then explore results, run bootstrap, etc.
load(working.data.loc)
load(working.data.loc.high.e)
# load(working.data.dir%&%'/working_fed_lands_data_modelrun_2021-11-01.RData')
# load(working.data.dir%&%'/working_fed_lands_data_modelrun_high_elasts_2021-11-01.RData')
# load(working.data.dir%&%'/working_fed_lands_data_modelrun_HighUS_ROW_elasts_2021-11-01.RData')
# load(working.data.dir%&%'/working_fed_lands_data_modelrun_high_elasts_2021-11-01.RData')
# load(working.data.dir%&%'/working_fed_lands_data_modelrun_high_prices_2021-11-01.RData')
# load(working.data.dir%&%'/working_fed_lands_data_modelrun_Domestic_gas_market_2021-11-06.RData')
# left off here
sum.table.num; deltas$sim.ban.lag.0
sim.ban.lag.0$sim$global.sim[,mean(US.gas.demand)]/sim.base$sim$global.sim[,mean(US.gas.demand)]-1

# Run simulations for the mesh grid of carbon adders, separate for oil and gas
source(script.dir%&%'/run_adder_mesh_grid.R')

# Run bootstrap for draws of drilling elasticity estimates
source(script.dir%&%'/bootstrap_sims.R')

# Run miscellaneous calculations to support additional results in body of paper, and model validation figure for appendix
source(script.dir%&%'/miscellaneous_calculations.R')

